{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "from sympy import *\n",
    "\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses SymPy to derive the equations of motion for a springy pendulum and a rigid pendulum.\n",
    "\n",
    "Here are the symbols we need:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y, vx, vy, ax, ay = symbols('x, y, vx, vy, ax, ay')\n",
    "length0, k, m, g, R, t = symbols('length0, k, m, g, R, t')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use vectors P, V, and A to represent position, velocity and acceleration.  The Vector class in modsim.py doesn't play nicely with SymPy, so I'll use NumPy arrays instead:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([x, y], dtype=object)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = np.array([x, y])\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([vx, vy], dtype=object)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "V = np.array([vx, vy])\n",
    "V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ax, ay], dtype=object)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([ax, ay])\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The only vector operations we need are `mag` and `hat`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mag(P):\n",
    "    return sqrt(P[0]**2 + P[1]**2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAAAaBAMAAADf4rbTAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAmSK7q0TNEFTdiWZ2\n7zJQnLHkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABsUlEQVQ4EYWSP0jDQBSHf00br03TKoi6CM3k\nagXBxUEcHMShg5MWKQpRB6GrLhYdFEEUBwcRrC4OCtZRXCoOBdFRXQQ7uRr8QwXFeEnvXiIt5KDN\n97736/EeDSBPyg46MgklTRgIkUxghAI9RMHQISJae1ez8GFfnjTLCTzCCEkPNCMiA0CoIBo7eDUE\n+h4xI/lG5aKke/Q22TKRUT5lAidE2PIm8ySiNVmwbUnAfh2Z4SlOsbIs9ZIkaG11jFLTrQ8pECPC\nucD/2bh7A3N6t/zDboezp1CLapUXgMwudUN5xASmAGXPaTjjLuGivIuX9dW8Yyhb3IRuJefWB4DJ\nBe7jRf51jLP8Cp5tmzM/4l6l0IlwOmHbFndbJb6kc1UeM05GHpFleEJrWcjENrBR5x+Zg2aa8w+m\n2caFWsNrVTR0/ofMupz8hbtovSF3C+VwJ6Jg74g6v8eabqEiLc2L8BD6yaaMFoMX6lfYUkpk5W4I\npdVvsr3FcYfZWOXmkqR3Lxu98l6xsDXtixDKefn7miOpfgwR+0AV84xgOePpg6rHjbSPQZ+89nEj\nZq+qjTLQ/AHXdX0I04TwoAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\sqrt{x^{2} + y^{2}}$$"
      ],
      "text/plain": [
       "   _________\n",
       "  ╱  2    2 \n",
       "╲╱  x  + y  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mag(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def hat(P):\n",
    "    return P / mag(P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([x/sqrt(x**2 + y**2), y/sqrt(x**2 + y**2)], dtype=object)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hat(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For convenience, I'll define intermediate variables like length:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFYAAAAaBAMAAADf4rbTAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAmSK7q0TNEFTdiWZ2\n7zJQnLHkAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABsUlEQVQ4EYWSP0jDQBSHf00br03TKoi6CM3k\nagXBxUEcHMShg5MWKQpRB6GrLhYdFEEUBwcRrC4OCtZRXCoOBdFRXQQ7uRr8QwXFeEnvXiIt5KDN\n97736/EeDSBPyg46MgklTRgIkUxghAI9RMHQISJae1ez8GFfnjTLCTzCCEkPNCMiA0CoIBo7eDUE\n+h4xI/lG5aKke/Q22TKRUT5lAidE2PIm8ySiNVmwbUnAfh2Z4SlOsbIs9ZIkaG11jFLTrQ8pECPC\nucD/2bh7A3N6t/zDboezp1CLapUXgMwudUN5xASmAGXPaTjjLuGivIuX9dW8Yyhb3IRuJefWB4DJ\nBe7jRf51jLP8Cp5tmzM/4l6l0IlwOmHbFndbJb6kc1UeM05GHpFleEJrWcjENrBR5x+Zg2aa8w+m\n2caFWsNrVTR0/ofMupz8hbtovSF3C+VwJ6Jg74g6v8eabqEiLc2L8BD6yaaMFoMX6lfYUkpk5W4I\npdVvsr3FcYfZWOXmkqR3Lxu98l6xsDXtixDKefn7miOpfgwR+0AV84xgOePpg6rHjbSPQZ+89nEj\nZq+qjTLQ/AHXdX0I04TwoAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\sqrt{x^{2} + y^{2}}$$"
      ],
      "text/plain": [
       "   _________\n",
       "  ╱  2    2 \n",
       "╲╱  x  + y  "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "length = mag(P)\n",
    "length"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`f_spring` is the force on the particle due to the spring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-k*x*(-length0 + sqrt(x**2 + y**2))/sqrt(x**2 + y**2),\n",
       "       -k*y*(-length0 + sqrt(x**2 + y**2))/sqrt(x**2 + y**2)], dtype=object)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_spring = -k * (length - length0) * hat(P)\n",
    "f_spring"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "xhat and yhat are unit vectors along the x and y axes: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "xhat = np.array([1, 0])\n",
    "yhat = np.array([0, 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I can write force due to gravity as a Vector"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, -g*m], dtype=object)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f_grav = -m * g * yhat\n",
    "f_grav"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To write $\\sum F = ma$, I'll define the left-hand side and right-hand sides of the equations separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "lhs = f_spring + f_grav"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = m * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I can make two equations, one for each component of F and A:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAA1BAMAAACZwnArAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMrvvmVREiWZ2\nIqtdv6urAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHBElEQVRoBe1YXYwTVRQ+s91uZ6btdgKRB4K0\nLv5GSFYwPqyLTvTFB83WB0KMIIs/iSLEimKQgFsTRWIkW9T4AESGNWYTWbVCIomSUCAKLhtoJBh/\naUOCRlxJSSAu6LKeO3PvnTvtTLd/a7qJN+n03O+c79x75t65c84AeLbrnvBUTTfFLdNtwp7zHfPU\nTDOFcnGaTdhzumrMU0UVwfRkFvXrpXz9PlqM4Es6zDnxw3DGw9lhD7yh8Nv1e4tu2/1eSjKC3b6Y\nuzNlpTveWDSUr9tf/4PoQtXUi3LW3dcs3R1vMHp73f76bkMXErTEvTz95qWoHA9MTNYKkNMq9+du\n2X0mRhSRDLm6tNYGnG2nXfwWQ+35YqTKvnRJGgsjp0cH95vSkqjSo4v57y5YMVT3HQsWpCsqyNpJ\nCLlHkssWj1l1PxyjlMGlrmNY8Pyq/ToJ7XG4pMKu5Ar4yalgvReYUPt/IGVx1Wx7zMULhT93UVUD\n+dOwPwNzR+bsSbvTVrvD1aAzqLEv2zruwqNwjsbrYtEQKHxJdHPE0eOaSfLp9dSwLRW+wjm2QOE2\nw4amQgoWRK/OnqWZmQUom09LndyF7H4QmrDrzuPM+oWAPQ10FoiVeuxHqCifvt5hpaZ515exRAmj\nF5oJhwoCMgWiLy46bcuIPUteC1CcT7/lsJpp9wapKGdsDCUTlt2eIYdZfZ1IQuTndLFnyV2Y7HQ6\nYWckx7kyyJw5I7Hg1qKF5awGCc53/0H0OntkE2xf/vOJDIA8dOyRc5dHoWX9geXieM5ISJhWrr0d\n3rXMWCSzN0J4JViw9LfowpTPDN2vWVTlk19GRpWhR0pMKgeiadH2HO6kxfAjJJZkJNwMH0IuIRcA\nogvgAsDwMDPlkWiIhDKY15m5duutWx62TFgkxmYIFRj8j6kb6CDNPEKUNbCCUlX/Moh+A5stek3X\nnCOS1QDfJaUXW5MLAY9nPHIiSXII9OuwBIKGnKdDsEhOGQj48Gfl2m0TExg2aTSScHIt+OMMvmrp\n7OusNKyj1Aew5sgZgHtixiLSbrStKpV6kmgpzSPshXorvk5+PXBck6QrgAdyWxx6sj4DAB/61eDT\ngZ2xNBJlby+SvyAOinJtGomENHv7/oWGjtanQTelatEU9GuAe4K3ydJrUU9IPeKahHoBrhEUtxQe\n/6g7AjgEdJEVimbhaaLD7XEX3R4y2fp3EFCYMKgdHfMe7+hIIIqHHqautJVE8hRI5E1sxopnPeZN\nuCdqbo7d5c9rYO4B3FKRtHYwC0thF5AVUnulHg3upsOw3QX4LU02TNCZa7PnBP18zadm7a4PyPIv\nugdR6SreMYmm6bjsXaCMY7fWFs0IzEjagFX4DANmFv2agfdpDAbgNbLRjM8wkm5qyyM5mYIWDcHi\nXJtF4s/DY9x/8dkljYE/tpVSVwHer0Dse25dtRAxBEpkWwYOg7QJMNsbCOrRbPs4vCJlcKNFUkZP\nli8+jyQahzcIvzjXZpEE4gqff7jkfTIf9uUNi4oJG96vFkOcDZ9Yn/hIuMrEtC3PCbgYG/C38UuA\ndwBOHwNl9KNOOD0CmE+rGyCqA8sPeSSBcVhA+MW5NotEOn8UD3OrtTI2A2Du8EOjKYuqdAJ+s5DP\n61xpC+G4LZeR/LEySvCTx5Y2X6r47AIIX1XyTC/+s0jETM4juxRp7nJ7ysLVdc+7GbDyTu1101qY\nbh5cTB9KyHkq72AYrGjJclkQlLTVOSt4KDeSQC0Vv6XQx3C2VAm8vCtzp9SLwuOKPgYHtRJPuy6U\nQCLQDct4F5/+2trLlLYY326lHnh555INMWv5j6NJJnv9+8gLwbvtPKpzZYSuEwcqFKQYNXyGHEQl\nzS7v7FOyxKgCIBSvwMgycbuhlZADSW7Vr3FRENi2qvdDgeCyvDhQXu2pFcof+j5zq+OQHk16+mis\n4tka3e3hPDVhifbJaPZZeReKccspFeTe2twLXwn2Uw/OSHh5B0/WNkK1LF+6WoZlH0KeNPLnziFQ\nDEU3MRaJWMcRxV5TO+WXWofx4cxmw77MQrhvy+vWE88icdRxaObPTHkUOEB4TfWjmInxCPJ2w6fa\nVnhuYsLyQSNx1nGoku6sfozqGTOTVXPCXYTSiT8NbhLZNBJnHUcMTmVFsymSH63e746bkRM0TKL1\nDYCIXnUc0YV1cp3aJqVq8N+fxq8EGiFiWS7UYOw5cdRxNfj/7yhtnQBvkuFexe+Xh+xxWSSOOs5W\nN6EUwormBpyXMuYvhHF9WGORiHUc0zXnv3QZ5AROTRo9NIzlH28sErGO48rmFPqyruUPi0Ss45oz\nAD6rnPE+lwXBrY4T1M0o+gvmVwKvqYl1nJdNk+DKtXy5mYh1XDm7ZtDdqzfDLBoxh68a4eR/H3Xe\ngX8Bc50nOQluB+sAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$- \\frac{k x}{\\sqrt{x^{2} + y^{2}}} \\left(- length_{0} + \\sqrt{x^{2} + y^{2}}\\right) = ax m$$"
      ],
      "text/plain": [
       "     ⎛              _________⎞        \n",
       "     ⎜             ╱  2    2 ⎟        \n",
       "-k⋅x⋅⎝-length₀ + ╲╱  x  + y  ⎠        \n",
       "─────────────────────────────── = ax⋅m\n",
       "             _________                \n",
       "            ╱  2    2                 \n",
       "          ╲╱  x  + y                  "
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq1 = Eq(lhs[0], rhs[0])\n",
    "eq1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAccAAAA1BAMAAADG/z0nAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMlSJdrsime9m\nq0Tz+RmlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHVElEQVRoBe1YbWxTVRh+7tre9q5d20AiBkJa\nh0pQYyYYFhMNjcYYf60/gKCG7ArG4UJio3EjE0L9BRIJV4kRo7g6RUmMbqgIJEQLGknwBws/nBrR\nKgZ/KNlIGDAyme/p/ei56znt1t0pNZ6kve993o9z3nPP13OAWsqNn9TiVW8+t9Vbg2tp72gtTnXm\nE71QZw2upblqsppXJF/NYuZ6ZXDmMSpEaNAjH6SUzgzeiEusnpHgnsK/ehptcrDEL4//UVADA9g6\nWWO9Rw9IFJ7CwVn9lO0PUWNfbUxiu6TRC1IShbfwHd6Gc0frv5Pe4+EMLrtx5+0HR6pZaJyoVkbQ\nK5stNdfKObbuTdLbcNx/kQM50e/B6ruHiycTm2ZxvCoXlVEN6ENwQFx9Q06MTwf9aQrGXnSmrJrI\niHJNBXYiJOnJXkPmOmVcS1qmnTviAid1/9sMXSxQeQQ1ZXGRktyP4bw44jtieDpoY8G0Vo2mpMDv\nKZxm6DcClUdQKI+TGWBhR7/kix2aeUVzrBA+wz8miHYP2ljdvQWBzmOoTxxP49ej4IqU0KoKk9ln\nOYUL2jWB/+dmfmFdoPMSakoqouqpisgIX88p/sWS59JnqMhklBbHKyBeqtvjZCEcyo6nB4IvE8qI\nwzSWWkgGos/dTribySiGK5aad159Vi2TLFqZQXDEMZsdIdr9riSwL8srRBOUEp/EZAJWKpbn3FKE\nTkt0W6g5BgdE89Wyn+VHLMdV4JqgNr4cUFvsl+LTnQI2O8pIzhLdFieLqN89HBwnmRDoeC577vvX\ntmTObSrIbKaIx/ivwiao0tFjRDe83j1EAVZ3r4l+d2UIDftOvMjFc6cA1gPztrzVlcE5/GaaWRYm\n+Ynq0RTByrip4/73djwQN20E2WxAb07fnlEuF4hdzKwk8pw/bao4a2h6Y2gbEgbUrHYVAUo8sQRn\ngK4uy7SUZJyQYIYS0COtvqT/llUvmCaWhUl+7l+1htlZh+d1zawUV7LoEeyCaVOeDa1hsZXppdAu\nTVoczRqm9d+b58zpGK8dxPzUg74kelOgRXEUxF/QniIOE9EDg6atk+RunQAf/dS4eiFghCcmRlwW\nJvn5YmKiiF4ydaX/BXk6pBRt/KVs5ixj5WaEs2h7mfYESpb1PZRFDF+aKrlPXWpLk60duC2FpqPd\nxxBPFEC7yWHQ5u7Ti4vuIfhSsDYIO8no8QFy/pp+Chqy9HCKZeEmP1cdtSX0x9FqEiSFy8ZStuXx\nLBtGtPxT37tLNc7D65knxSqV54HYIHulbeM0/OMItoDyxXLQkkTj9yCp1ObmRR83N+dIRIBNs7uY\nBG5ucxYu8lOW5GEodBQp2pRnc8rADjaMfDqGjTirovbiGq47FcR0Fou2jXElMoZQDsNg3awOKG1x\n3GvWY39JgO50A0UH0CDgWmJbuMiPOVzXs1G37D6KpFyiL6WYBKk8G+rnUdApiTq5T0ubFdf6n8hw\nnuMqC4sMPqOpoFLvxtJYh5W06Ib1Y5RkcVOnvByfrQU0sNwC8a0IMsEqtoWL/ExeXZVRhJKPmASp\nPJuE0TTGzoM0bW4lhjGjYn45K8S1eQgmMd+gI6ia/EoZw88GPlQydCCLFfQ2A9ZZwU6Blt0sHmW+\nw+ldeJNrh23Bkx+tbJ9cjC8HdZMglWcTHXqyBWeBLcDaJ7jQbpGffUJ5hNmHBzmvHh3o2ZxHtAXa\ntjTOd60H9nSDmIz6HhIpWOdfOwWgcQxLmP/C7nkb80ywim3Bkx9/2el5YdfDQwWTIImyoblSrWjZ\nahZFfShZyWwHp/QVJq+uoE0syneSY20nyZ+GJYd38hGdmJEqLnlORLHQVDBxMVuHDasDYndC/XBR\nl2DO3iejecdnV4PhyJxgWbjIj6Qml00phHqBJmvV8oppIWHrDizvXzwdd1+ZdHbGy2odPlMGcYCL\n/LiDOVYuGwdF4MeedOlNJn1kKiRs3YEFJ0o74toTNhu2kfKnj+fcZWoX+Ynly/QMcNkILeSgkjR1\nErZegqcwKOS1IJitoHSraHX2ujQ6H1syIG14Fu+X3Dmtc7968VZishK2bsMJpze8qLZCjKMVdDWq\nNjp+nZZUWtaLgA3T7v+PlBmTwvJWKi02FslZkjtJB8antuXsPn15z+MHKaSYreOG96EdKJF4HPe8\ncmHAWajGJ2Xr0B+ji7MSiUcoI2yUx6B2xNOACovWLWXrWroPoWyJxEO529PqJcHmerq+aXS7Bna5\nJGHrCnEJjuGS7W5D0jAv4Ze8DIY/b6dwEZ3F5HLh2DpdoralmNouWsqWZu+pFLyN3Z6ny6U4i0m5\nFJ9mfHt1JRJO1xx1XsI0VFdTDjK2TgflmR3lrof+CRI5vYkaImPrjdno5IuI66HZ02uDcgWBHLnI\n2Lrybc8Yqeu89BtCJmvPSbq0SNZ5htT8Xv28KAmbz5+eyr2CKMB1hYVGllRqTyu2VVLXiS7612Cl\nlv7ek6qkrhfdiv9EFlV6e1MV/f/qf7kH/gZBD0jIsHngkAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- g m - \\frac{k y}{\\sqrt{x^{2} + y^{2}}} \\left(- length_{0} + \\sqrt{x^{2} + y^{2}}\\right) = ay m$$"
      ],
      "text/plain": [
       "           ⎛              _________⎞       \n",
       "           ⎜             ╱  2    2 ⎟       \n",
       "       k⋅y⋅⎝-length₀ + ╲╱  x  + y  ⎠       \n",
       "-g⋅m - ───────────────────────────── = ay⋅m\n",
       "                   _________               \n",
       "                  ╱  2    2                \n",
       "                ╲╱  x  + y                 "
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq2 = Eq(lhs[1], rhs[1])\n",
    "eq2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I want equations that are explicit in ax and ay.  In this case I can get them easily by dividing through by m.  But for the rigid pendulum we will need to use solve, so I want to demonstrate that here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "soln = solve([eq1, eq2], [ax, ay])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can extract the expressions for ax and ay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKAAAAA1BAMAAAA0Sv3uAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAELvv3c2ZVESJZjJ2\nIqu2f7MxAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEpElEQVRYCe1XXWgcZRQ9s7vZmdnszyAKBmN3\nrVCk1LQPPjQouKAFxQSjkmoRYa0GlLZk1YeiaLOiL6nSTNUErFH3QQKKkH0RX4RMU02trOmiSFXU\nLBUflBQ2QknqJo3nm5/d2ZkhZu1DH/TCztx7zv3ufn8z3xmAJt2c5TV+f0kEXrNIL9oaX3Nba4yO\nnACmPShDqR826efcyF53QD+VF8CQB2Uo1xzSz7mRZXdAf0wTwB0elGF0wCH9nAtRllyBcMfFJVEX\n11aL6MBbrVBQpPJ/W2wfpm5HrAZ0lQ9h9Naf5g2OdvJ05tjOHgN7X5itWNlS+cLrky0N7aBDjz1R\nQvf892fYkCatdF2/G8kMlH78gNxWQ1oD3sBCDiMkL1aSVcQn2dcufGDsMVt4LumXj79SkPRYX8Tq\nqlw3mBE28G1eOtCZ3yNGLy+JpeI6yXWEMvgQR4DjeF97yVPLDEc+403V1CXZGkv0YoHATAm/fvyV\nJkkrYvThDGYquNdcF/7VfUhp0PCg2d53GX6EkISOjM1Evv6N3ingsgC4VZIDmCnic84FwHVZKHEO\nwnlyqyLBb32/mENNGTaVLqXzFfRIuCQAbpVUUZuuYBuUGrR0AW/ihiVEikBnnf3wm1SXlhOEOUbN\nZKeRLuj4S0UvH0CEdYxoOpdjGeqAXOLz0xsTBQ08E6/hpL8eZ0haUSFrZxG3Cp5HqqRjpRtfQDqE\nhQLGY6V0JbmGaLUbnI0VtZNDNpTlUC1RDCjI/VFXMZbfiR8tdjtiB4E5HfGDHwGvAt+chrL4dhaJ\nHQa2AxO6VEeqIC2ePEPeb6EipgxsKXefKPrJJhLKNX3cI1b5iqzEeXUVmBP78EpMXUKPu7369OPu\nsH1f/mMu336r/1tc5Rm47m5hj/p6sf5vzFPFXcJD/VfDYfecBPrtzUwi017+P2YnC1aKuj/wfTCx\nrd332Hf2X74rXtY+Uys8wdqzp+z0fnGO+ixS4XHQlklOB3aL08Vn4UKCZ2s7Fs03ske0hutyqCna\nsmub2X22e2MTosfz1GejFyiwRmcLPoLAiQao5mz3xQYknImWyAp0CqzVgrwrgJKyDXDK8VoKxnIO\n3LybAuuSKQihvPdzeVGZvMVh40XKHlMHKrpSslC7oDRBfaaN4jUnt3EXAotTywOfIiy0A+kv8ZxD\nRuhYOvCTw8/ai2IXVDmks50PH77JyW3eKbCiWSEIgU+pERd08wtAEgll/iwdeOf6ugBodsFzlFFb\nw+vrNQt1X8lQpY1V2AOhrrg9zlNjUwwCWf68OtAuqLEHqyLHbxRYrDOUyJOizHrS/AI4+hCjmG5m\nu9qNDw7eNThofoiMadRzgcZHgCptnyrIIVCnKmsc70iRm4yd9ulAu4dMje8StN8osOaBI+8Iplfo\n1OjAOYrsLPC8gLw60CnYg1BV8Bsan0wq6g6dQ43zIX2AyT4d6BTcj7HihsUEqWQRr0L+vcTt9yfk\nHCGfDnQKbpkcrog2m7bhSkdgA6dg4OfgRtUX9GOB9FELTQ6Ir4J2LFR7bKP0iBEyNuL9nHK56geb\niFI+0Aw25wV/02+ubWDWbCB6dcG/AVt5k8Zzd9rPAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{k length_{0} x}{m \\sqrt{x^{2} + y^{2}}} - \\frac{k x}{m}$$"
      ],
      "text/plain": [
       " k⋅length₀⋅x     k⋅x\n",
       "────────────── - ───\n",
       "     _________    m \n",
       "    ╱  2    2       \n",
       "m⋅╲╱  x  + y        "
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "soln[ax]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAANUAAAA1BAMAAAAkKrMyAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMlSJdrsime9m\nq0Tz+RmlAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAE8klEQVRYCb1YXWgcVRT+bjY7+5vZRSX1yV1T\nFa0gIUpDXnQQxMesDxWlSqYVGo2KW6Wmxh/Wp2htdBD8QyRrFWMRDP4Q39oFwUIoJPjQVPRhRWkf\ntLgpGzXQJp47sztzZ+7cndgmOZDZc873nXvm/syZOQE6C3ujnwiJe4wwmgOGIZ7v+m88PUrrqnLG\ngkxjQ4ADypjo2SUanfVciePHZFK8ATigjImeFdHorP+V5/h3MildBBxQxgSPviwYEep+jqeaMitm\nAo/J7qBHKwY9avsWzH+LLK1X7/hxnHv9vZfLQHz0UOXCW1Nl7HrksOWE9s4g9VXYKF1m9gtDifoi\n2FrvDSfQU4E+hHdRfa3MVoHHMV3FCMAuWj2LSIyagLkXCbohWQq/7vujrkR9/HiT5oFMGR+W2Kfd\npQG+nPFlfijouMSbSFbwAx5FqnSMqyEycj9tgRL1BaQv1skeNvDzqZfyjK3x5cxUMGxhN0DHg+7i\nLuTyjE5Pjt+VJLN30PyVqI8ee+FHsl8ELnM3nfOeIoZreA6UFjET0wYtaqYEOm50QyEy+FERatQX\nUDAKJQtTDH9zN00kV8svWDgCvYF8oY4ncd0yYjWO0A3JwppsJaVE/fwFFOomLmmYo0qFjImRvEmn\nYgVaMW5QNZnL8lxlJBcx5Y90rGyDrWlK1B9xFjnDxNoOPAt2HNN17M8aBatnFenFHaDlXdO6aQ3L\nSFf0S/5Ix6IT3NSUqD/iKLIzwISJxMxp4Dfg/UPQlw70IzVZxlFgzGRN5OpgZyboYZAlWcN8WYnK\nfNmTrAq+u+kckslrllo6o+o4g/ZQQCfo+aLl9PkEmKtKlFkBZsDUlv3HQPvsY2AQkwGaaCrReOgz\n6YXGf5koeVZLOz9hSD7PoUSjcnlDXLHGi7cdvA257OLNk219Lqd4b22ua+7kclOreG9KLraTDzlg\n8MFCxS7egNbXt/Prvr4q56xfiQQGF4ZouJBdvG1r6/fLLt7blMsu3tuUyy7ekblmhfUPV+0hIi/t\n4t1hv1KVyFE2QvCKt15T8nvqDqR9TuVWlrEjedkpe4LFW2Zwz9st9xP8lSyJZtGHzwYktHhLcV+2\nPEP8o02SmEXfA5slrNga6QT/8pAkU0/Rh9wmSbrkDjQSujNUVDdLrvUGGnTUwDucPt6Coj/9wfiS\nPvpm0B9lP+UStKqjBp6PMZfgKlpyEoVXsNd1bExh/S5vvqX5c2WrLsFV7osVMW2GNcEuhSvB6Sdq\n1AaM/3R+FLqpG5zhvledBuscfre94oV/iNPmnhV9IXpw+jHi9OL78gDu3fOQczba87IbrO6b97wq\nD0NF/ZPQJthHFabPODBOf/vwTP5hnFxfd5itXE6DlVlfbzhu8Up92m7oq/YIot+ve9NPUbMF8O3K\n41aR1MrVqcGa4+1TuviOGBaiu9P/83ZCs6ZN+cclCu9wZfsFeryzDXSZTrAbKyne9Edq1NnZW9Td\npDbRlfZ+qdovOmD9SCwifsZwY8IVb/oZWr8HOemBRAMHPXY7l6r98pgRmjD9BBW6G4muryQbqZoX\n186lar88ZoQmTJ/9i3iV6Gzp4POn6bct7Vyq9qvN+1+/s1aXFRLQzhXVfoWEql3T5oUwsP0OVzZY\nYUFRvmTjtk4UZYPVKUiF6ZcXVRD3KxusTkFKLPw/o0r6VQGHryp624P/A2znpUtAy0gUAAAAAElF\nTkSuQmCC\n",
      "text/latex": [
       "$$- g + \\frac{k length_{0} y}{m \\sqrt{x^{2} + y^{2}}} - \\frac{k y}{m}$$"
      ],
      "text/plain": [
       "      k⋅length₀⋅y     k⋅y\n",
       "-g + ────────────── - ───\n",
       "          _________    m \n",
       "         ╱  2    2       \n",
       "     m⋅╲╱  x  + y        "
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "soln[ay]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can get SymPy to format the result for LaTeX:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{k length_{0} x}{m \\sqrt{x^{2} + y^{2}}} - \\frac{k x}{m}\n"
     ]
    }
   ],
   "source": [
    "print(latex(soln[ax]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- g + \\frac{k length_{0} y}{m \\sqrt{x^{2} + y^{2}}} - \\frac{k y}{m}\n"
     ]
    }
   ],
   "source": [
    "print(latex(soln[ay]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or generate Python code we can paste into a slope function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k = Symbol('k')\n",
      "length0 = Symbol('length0')\n",
      "x = Symbol('x')\n",
      "m = Symbol('m')\n",
      "y = Symbol('y')\n",
      "e = k*length0*x/(m*sqrt(x**2 + y**2)) - k*x/m\n"
     ]
    }
   ],
   "source": [
    "print(python(soln[ax]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g = Symbol('g')\n",
      "k = Symbol('k')\n",
      "length0 = Symbol('length0')\n",
      "y = Symbol('y')\n",
      "m = Symbol('m')\n",
      "x = Symbol('x')\n",
      "e = -g + k*length0*y/(m*sqrt(x**2 + y**2)) - k*y/m\n"
     ]
    }
   ],
   "source": [
    "print(python(soln[ay]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see these equations run, see pendulum.ipynb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rigid pendulum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving the rigid pendulum is almost the same, except we need a third equation to represent the geometric constraint.  The simplest form of the constraint is:\n",
    "\n",
    "$ x^2 + y ^2 = length$\n",
    "\n",
    "But this equation doesn't involve vx, vy, ax, and ay, so it's not much help.  However, if we take the time derivative of both sides, we get\n",
    "\n",
    "$ 2 x vx + 2 y vy = 0 $\n",
    "\n",
    "And if we take the time derivative one more time (and divide through by 2), we get\n",
    "\n",
    "$ x ax + y ay + vx^2 + vy^2 = 0 $\n",
    "\n",
    "And that's just what we need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAAYBAMAAADQRaYKAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIquJdjLdEETvu2aZ\nVM0GsGrEAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADDUlEQVRIDa2Wy2sTURTGv+mkuZM0qcGdi2Kl\nXZcRqi4sddBKFQrmL7APiqKCRAtuU7AuRKRdKAqCliIIItKVIAoW1Apd1BEVUSjNyoXoEF9UUInn\nvjKPZKK2mcWdc87vfN/JnUenQHMOc+lZvFHfnUI83DjZjp5YEzPXUoqFTQCjGM7F2SRzrZ/jWDPq\nV7F1Ms4nNWn8iGPNqa82uKfZtebMiHMZjwNUTzoN4MaRudDAoy/K2OOhfLS2/vxhA2nbQhS2u2w0\nWqvmLPappZY60LIttZFTVY9qcBbnqrEM7gP9kZKfZh0/ronqwMMr79Sz9qqmvfXTyq1IcQ6Ifzbr\n2PvyOvBGpaJ47exUpVL2xSL6Aky7gNG1/yjYwLapXrUKqO0FnXgC4wrrc7BF7k1Dqp+Zl1LfW82m\ndnYCQukzHbFfNNsBzufSZUzgnnNSraJB2wtqv0WibGZncEiKNaR6slNKtSugZvP2MQilZv0ePz5S\nanwFhm2wa8jMYhf2FF6oVbQqe0ENdxnp2Y5MCbeljYKGu4isLaXavzq7IzOCl1LpMx2p2Zk1/uYX\n8IHqchUNyl5QhkvY5BRSDr5Jsf5hVLfcgEhQtW/eflEqpSi4Mto3XfPULIp5qktbuZqe9/665y1A\nUWsNw3lMF1q/U6MPQfVEVUoBvTied0xeVmqHwztIWXvQ/V51UbSxShsmWwa18la1NUnpah8ElpGY\nkS4KIjMCMygSVO0bi4ANqZQq7DzCjzGeHQceFVB0cTebe5MoYxByFZ16tqDpTtALOQc6ByHPXwdF\nAurZc8jmeQcpaw+xFbQ7bWOm9TtdNublKhvVbEHpYbR+AkuYng9BmCVmB0WC6tn99NooZe3spMsu\n09/H7h2bD7DuwX0PKORryF5QsK4h+vif7n2eC0GgKywSVM9u2U2/VCqVZ+DEevfmA2k41LdUVTMl\nHiyrLAJVVZ30bJVKZbjlL5mlLq9o60FxEi0lpv/7CMGo0YVgQSiDhf+Ox3GTf//TzvqU/60KCqb4\nh94aeBqs/VsslHGtfwB1oelRKPxU1AAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$ax x + ay y + vx^{2} + vy^{2} = 0$$"
      ],
      "text/plain": [
       "                2     2    \n",
       "ax⋅x + ay⋅y + vx  + vy  = 0"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq3 = Eq(x*ax + y*ay + vx**2 + vy**2, 0)\n",
    "eq3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can represent the force due to tension as a vector with unknown magnitude, R, and direction opposite P. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f_tension = -R * hat(P)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can write $\\sum F = ma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "lhs = f_grav + f_tension"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rhs = m * A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And make one equation for each dimension:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAA0CAMAAAAdSEKiAAAAM1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAADxgEwMAAAAEHRSTlMAEM3dMmaruyJ2\n71SZiUSDyi93ngAAAAlwSFlzAAAOxAAADsQBlSsOGwAAA3RJREFUaAXtWuuaoyAMBVQEuQzv/7Sb\nBC/AUJWOXd39hh+tSkJOY4BzmGHsjTbIECS0cXjD+R6XMGJcFb/ugdAUtQ+a7GXom/xuM9Yz0C6o\n2zA0BR4nMjfBNrndZxwkxjYdFgbXWiqvXSyS+zDtRe6D1HoMtEZwx9gAyDux53Fznw4cEAwTFrCD\nSwf35snVPJcwVTBCF7TI3ZzFvfCxhFX8AsN5kdtzubWvj9W7rhE9LXKY6oe2eRV2mFhuONMBgOoH\nAxYIECYdAO6/gmEj3Kvn0gorgPkQPDtqzWGB49oNz8VLyf39+M3Abwb+hwwo6XTf4XKeEEDvpFLA\nBzXzTo/+lt/pRNI2Qu2RS+GSzlIC6Ng44v4UNIA1kYbfgvpbUD7hUm6RRiUEsDeRt6oALBa2Kvjg\nolubIIUWPtQg2E6zRFmjJACzmQByyC3WgaFe+SRWOCEYvoralQAa4gG2w9863aNq3PpG4YK0FkHF\nd47oiERtBNCSkumw14PkfQzF4kCiYM4JVC4ZAZwQaqwLwK7Q6ky7oKwPwgh43QMk3EEiEwKoqIRp\ntjHhoPdcU5+vdqB9rlcSVq+MAHoqXh1VunRnz23kPSv2uWxWrFQh4Z0d7fnqb7OuhG9/ZPMEO9iS\nZPEbXg/aZv16nIYeXoDrAHB/+lStzboB1mtTW6wlEyQ8zt7XPltPm/Xm9/4Vp4la+Me9snj48rbN\n+uUwJztcTWDSpp8MwHdPqkrrxPH6y2qC3UYLY0Szt81/s74eZTKiriR4KPGyPcDfrZPhr7+cKzhR\nAUD+gYnkm04BmDuLKVcTLNcV6+tRbqxIx3edqoBeGGPGfOcoALuoDvAQq2Z9OWBFdI6GnROcqAA2\nERXKg+aAgVQRLRyBg9Ssc9+f3/lBLOppmPkR5nPvGDgHDNZEtlca/nNMByOYeKLKGKnY2bgePmrK\nbiJluU5Gj5om0vCDUBd1T7EUTEJANxVQiZFnGNQkbud6eU8Vh6sfaSLOTCxTK1MBlWAl4BEZLH1U\njD/xiNNfZc3yijMVUItXAka97pe6qjlc/swiHRPLhpupgFqsEjCHYxwaomb8kWd4dOEXTXscoQSM\nHlHwHvteZAGnRmuCj4f0+QaO/D5K3WPXqyw8/EvE22MhPcMzsb/auvUYpj2s19rmZKN9jGYP/Xlt\n34zpH3T4A9v4IFjKd628AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- \\frac{R x}{\\sqrt{x^{2} + y^{2}}} = ax m$$"
      ],
      "text/plain": [
       "   -R⋅x            \n",
       "──────────── = ax⋅m\n",
       "   _________       \n",
       "  ╱  2    2        \n",
       "╲╱  x  + y         "
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq4 = Eq(lhs[0], rhs[0])\n",
    "eq4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOYAAAA0BAMAAABsslmsAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMmaruyJ271SZ\niUSa32RUAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADxElEQVRYCdVXS2gTURQ906SZ/GY6KIgbSQyC\nooJVQaxQDG5dNCC6EKGl4CeLQhRUKiKDCxG7yBSFUq1aFy5cqHWltkhjFYobG8VfodIurFBBTEFr\npWC882kyk85LIp0M9EIy93PuOZm89+YmANMO5Y8dO82s1qjwG0inasTNoPX+AXw/GcUapQOk55mv\nETmDlu8HhjOMYo3SDW0INYGbUHBFqpHEEtqWZM8vQAj14/aSWq0StyUcltEVjuJsrSSW8J4DGjKQ\nggpoA7tkc3Q824Afkte1A+OlxVS37QgC/S7dJdTjOalwOAOe7tYdq8sBrcp+jOJHwh1FhLfNJRDu\nVLBmfEx2SdMkM2LyXXF9Ue6vK0ImEY/CK6bQFVfsfemKzkoQCY3fSIkfLvcmxfELpZ/X33FtO9a9\ngv8OC1HaUV38EemMwHcicguDaocQ0yyj+t9kPofsIAI5M0KtLMtCs2gY2O+JIp1FawkTdxfhlH9g\nBHyqiFi1U7UNQP4/TCUuNCKYQossRaYwLOF4iWZ4Fh56gN5Dg8JAlDRUGbYkcJXmAfACRG41+jyR\nOMRZtMQZCCu+2qhVxnl1/qAZ4jxHXYL23e1sIzeSJSmJ5v51C6JaaiaO7nABuA96YIWjXVZYZACn\nQjJNpjMshI4fq7iyJbyybx7+v6jPoS6btdbqlPpdAm0j8RcLoeH9KWtbxUhMfm6E2IgA/VDsiFvh\nXPLdpVvgOvrmWQgN75uytlUR8ZlKIFrRcnbRKAqjz+1gE+elknQcdE7KWnslxGuj/RPabYgE2Re1\npoVZdYOUtSZ0lq1zUaO8m0664ZouHtlb8h8mdLRvwASwc2f64nbpQi68SPAA6alCtuAEp2iHOm2r\ni4T0LLMxerw6be+LhE2Gu7aYIs+jWEIHAq6xQCJkDPdpIaU6E5bIiSCQKLAcWfQsmvWZxbRjVw8x\ncb0nZsYhZsW4TmtoatMe3/HFMTGSUrl66bUOB5Qd2Df02NhDhqY27b1bhk46p+lvVrnU5XyDSekJ\nHubzBrmuqU/7YD6fM9IOXKY3E0l9lt4kbLTw6Zr6tLcUlh8MJ+ggSBrPnyLb21hsTyy2lRL6tC9W\nnPCCjcAzjYj+5WqLa7Aa66lPeyeUihwBeqitV8NHgRy6i3kYmvq0N+UdcLnfCGWIR1zgc376ngu2\neJ/qtHfaxuQ6mTi5ZPfNg2ZuQ1Ob9ua8E346+9WWxtAE/cSyrS8nyec22bZP69mK0962uUJSnGsr\nh6g07cv1smt74+waUGnal+tl13rYpRVe+Qf1aD67A/KKzAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- \\frac{R y}{\\sqrt{x^{2} + y^{2}}} - g m = ay m$$"
      ],
      "text/plain": [
       "      R⋅y                  \n",
       "- ──────────── - g⋅m = ay⋅m\n",
       "     _________             \n",
       "    ╱  2    2              \n",
       "  ╲╱  x  + y               "
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eq5 = Eq(lhs[1], rhs[1])\n",
    "eq5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have three equations in three unknowns:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "soln = solve([eq3, eq4, eq5], [ax, ay, R])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can get explicit expressions for ax and ay"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOEAAAApBAMAAAAxJlOQAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsy\nme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADmklEQVRYCbVWPWhTURT+XpL3mr4maWih2KmP\nCiJCbXGoqMtzEKGDLQ4OIjQUsYgUMumklipUByHo0ohDBnGoYqMdpZjFwSm1U0sFq+DgpBFFECGe\nm3fv+8u97WvSXui753zn57vvJ+crEFw9B88csoLQ/nparmMuWdhfjmB3I2vUTDuI7a+nobPcBsPM\nSlZZbWx+kMe6LTkeBTXsuPqNXMYJeY9qBdJzvpOnB9CknfoVAPzOXVTtmB9wbDM7jJiMMVVuTm5C\nMhOJ302gAL5ifEIfEZ67T5be4qzr+YxO2TF8cW6atWbMRZaymHYdYfQO9kznhePf1/yO2k5a6hjm\ngGS0g1OTxHZn95HM+OywaawCMfqLtlh2hNWxXdpNapD4EaFLIyVTipR5AZeUeXpOrwDzyngoUM2G\nAKmberz4Shpg4PvF59RkUxkPBW6FfLmbqdfVT+17vU5VR+WVzegGg8ypI+W+NSQ+B+LajKXNQ6E5\nlH4gTxk4xx/SZKSHRQRvGMk1jK/mXiIWvA/D3MIdleZQenKEZQzzU9JIi7YeURr9urufljbQVQ7U\njKVn8UyhOYnSOszcWLrgHJnqxidQ32k17ugBJWfKqJ7GFzRmff89tp4QnM1YuBjQHG2Bxe5XAI3S\n9RLL+EmZbI2XGtvOF3aP1TwOQ68h/GDo1VhwztHUiNJpek9mU/94iO4x2mKMt22sgB5Q+HNbB3J0\nHqnm0BM3gA3EtjhP1eYGFFoq4FOUtwT8RdcIlkUN35dhVlSaQ+kvQBW0O8v9VhVa6sLsvgbs+C+k\ny/ofXiy2q+iDSnOMgpZjP/vJPM92f9gKLXVhdjT92JVZaKNDYdGNk9ooNWf0BlH1Tn2yOaM7cxRa\n6sKZxgvvWqU6epO7X40JwsoeerUKLeVwOgdUMDBB/6vQZZcrXtDEvwaaT/WSltNHswP9OJzaglFj\n38wcXgfiUZyk1WXxvI4Rr2CGm6YIOr6AT8I8PlQCzg9VnMAurvrgR5Gd9KpdLQ0yunC/qGlvv87K\nnanvaqlgdETCg9lH0/ZKlKmFM/U9LRWMDZHwYCy2TUcN0lm6OFPf01LOmGiIhAfvBR/vEZj6hHFG\nNvW7rT3k8bfyNTaKxYVvxSJ7ZxKR8Be1ZQenvniPEpFoi8UrDk99wSgRCa+oLSs89QWjRCTa4vGK\nw1NfMEpEwivaU0swtigSLZxFzztFrYhEC3S+klZEwlfegtmKSESi+Q9/BPjo4L7BkwAAAABJRU5E\nrkJggg==\n",
      "text/latex": [
       "$$\\frac{x}{x^{2} + y^{2}} \\left(g y - vx^{2} - vy^{2}\\right)$$"
      ],
      "text/plain": [
       "  ⎛        2     2⎞\n",
       "x⋅⎝g⋅y - vx  - vy ⎠\n",
       "───────────────────\n",
       "       2    2      \n",
       "      x  + y       "
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "soln[ax]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAR4AAAAuBAMAAAAM6mv8AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMiK7mat272ZE\niVTFV83YAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAExUlEQVRYCe1XTWhcVRQ+LzNv/t/royK0i9Jx\nFHUl9Q8b3QQXLkTIrFyoNOOmiaumPxS7UENBwUXt6ypGCEldCEKxAy7828zCRYtoxoqE4sIRtP5A\nYwqGIpWM55z77n333nlvXjadZJELc9455/vu3C/33bnnBCBtjAVpyLbk3WM7So8zfX5H6QFY2tUz\n9GDu7s/Q7dk9P8O3Z2v7450+nvE1Ebw6nf5zHYbF355ynvOdmALwI5zVQ91/QQ+8Zqmux7rPWEgZ\nNjqk+yl6ntM58CRMNo2ECvIt5aKTa+Y39Fj3GbtCGTY6pPspeu7TOXAJpjpGQgVm+at03P8UZDmM\neW3MsrFQGb5++TBR7FG09/TvwKaI+IyV9m9ZCS0k7FGK2WhAtjtlLz+ePMe1l8/Zf4g2j7B5itlo\nQLa7alG8WSsRhQN5e6I+jbApSrDRkUz/E4vxWRQ7TROo9My4KnXbRKQxlmujx8acODzKW6+h1q1N\n8Aw/NCdOBmb8LXwtEjYRs4wVW+ixEbyt2eqyyXt27YZY2F7mS5OXf2TtVZGxiQACK9M3szFnDo+K\ndcJvLtwAZ+Htq9fhcr8vJshlDpwB9yKmTuAH/XvasP/9VxbDSr+/bhCd1dB5hhjuRYHVaOfZCN7W\nbKGFPK/l3oYD8Gn4WDxJ6un+BWVa+Sh+0M/NOd3qeK4+QPT8ZThMDGYjnP9Hmpib7Y3NIQevnE34\nHn4OfosnRHrc3glgzQ8DuL2T4He9wLvlNweIR4qH4KZiI8z3JZq9j9N4AKCfNejPzs2i+Qjw8g/g\nwXgRgEiPAx/DnhCBpwDIr/UcGGslEINKCN8pNhLcTWl0dpZP+5O/A+VDSPxXkb1G4/4PGw3Siidg\ncgKfuD/kl/HB+vCJQyOuBFg+JRshuT/M27IhPdUNKODS+MKdeJ48P8U6fEBZ0oOvxMMH6gsoJYYk\nnsTzBZKNENfbhKKb3IzIPqfUwlezDnt68AsexAvRGviQyxTmYIay7+AH/V/BD5agnKBnBvwJYjAb\nyfzTYkPT1UhpVGSfQxeEswHnm7XNwrrbVtOUnmKrdoeytEle3enCSu8YvBzzFPEa/kLx/hNsxPlm\ns6+31EZF9jn8B1xZ/AGcdy8sfpWwjHNunhsdPB4A55Cxb2H/T5pupadEWclGbsr9nNKoqD7nEmuY\n1pQIV74v/OI6ZSodsgkjJjIo2OiOhdJYk1IaFbx0aLyBJw+cwd6qFu3BWTjYIV6xSzZhSCJDko3B\nQfp+NtYk2ajca+ajPud5wHsQj2HaGIfXGMovpzH0vGRj7nPKs9EJ6MtG5U8j782KsNSDP75404CM\n4Or8hIjfM9IpgWIDPEEUNiZXNSqmHtXn1E16arQ3FUkEqvR+2SiYC7FoRigX6cFCDC8Fqs+B3xV/\nuONHGzqcpdAXyWMjU6IQx41KpIcK8RKoPgf8npyQ8VzLwE34OIVsZF4U4rhRifQcwV/k0bjPkey7\n/7QLcaSHCnFcOe++DG0FrRBfazSebjQeInAl4FZJ443KNQtxtD/YeVLRGv2wC7HUM0NVeBuGXYil\nntOw0t4GOQOFWOrZd/2t5nbosdeUesT/LjY6+vgbsWSpnlDNR69GrZgLC6EKdoBTWzg1ShX/AzW7\nVy7awE/eAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- \\frac{1}{x^{2} + y^{2}} \\left(g x^{2} + y \\left(vx^{2} + vy^{2}\\right)\\right)$$"
      ],
      "text/plain": [
       " ⎛   2     ⎛  2     2⎞⎞ \n",
       "-⎝g⋅x  + y⋅⎝vx  + vy ⎠⎠ \n",
       "────────────────────────\n",
       "         2    2         \n",
       "        x  + y          "
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "soln[ay]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can get the results in LaTeX or Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\\frac{x}{x^{2} + y^{2}} \\left(g y - vx^{2} - vy^{2}\\right)\n"
     ]
    }
   ],
   "source": [
    "print(latex(soln[ax]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- \\frac{1}{x^{2} + y^{2}} \\left(g x^{2} + y \\left(vx^{2} + vy^{2}\\right)\\right)\n"
     ]
    }
   ],
   "source": [
    "print(latex(soln[ay]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see these equations run, see pendulum2.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x = Symbol('x')\n",
      "g = Symbol('g')\n",
      "y = Symbol('y')\n",
      "vx = Symbol('vx')\n",
      "vy = Symbol('vy')\n",
      "e = x*(g*y - vx**2 - vy**2)/(x**2 + y**2)\n"
     ]
    }
   ],
   "source": [
    "print(python(soln[ax]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "g = Symbol('g')\n",
      "x = Symbol('x')\n",
      "y = Symbol('y')\n",
      "vx = Symbol('vx')\n",
      "vy = Symbol('vy')\n",
      "e = -(g*x**2 + y*(vx**2 + vy**2))/(x**2 + y**2)\n"
     ]
    }
   ],
   "source": [
    "print(python(soln[ay]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
